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We study the probability density function for the fluctuations of the magnetic order parameter in 
the low temperature phase of the XY model of finite size. In two-dimensions this system is critical 
over the whole of the low temperature phase. It is shown analytically and without recourse to the 
scaling hypothesis that, in this case, the distribution is non-Gaussian and of universal form, inde- 
pendent of both system size and critical exponent n. An exact expression for the generating function 
of the distribution is obtained, which is transformed and compared with numerical data from high 
resolution molecular dynamics and Monte Carlo simulations. The asymptotes of the distribution are 
calculated and found to be of exponential and double exponential form. The calculated distribution 
is fitted to three standard functions: a generalisation of Gumbel's first asymptote distribution from 
the theory of extremal statistics, a generalised log-normal distribution, and a y 2 distribution. The 
calculation is extended to general dimension and an exponential tail is found in all dimensions less 
than four, despite the fact that critical fluctuations are limited to D — 2. These results are discussed 
in the light of similar behaviour observed in models of interface growth and for dissipative systems 
driven into a non-equilibrium steady state. 
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I. INTRODUCTION 



A. Motivation for the Present Work 



The fluctuations in a global measure of a many body system are often assumed to be of Gaussian form about 
the mean value jjj . This assumption is nearly always true if the system in question can be divided into statistically 
independent microscopic or mesoscopic elements 0, as dictated by the central limit theorem (see App. |A|). However, 
in correlated systems, where this is not the case, there is no universal reason to expect the central limit theorem 
to apply. The fluctuations can then take on a multitude of different mathematical forms, including those of other, 
well-defined, limit distributions. 

In this context, the most studied correlated systems are critical systems. At the critical point of a second order 
phase transition, a correlation length, £, diverges from the microscopic scale (taken as unity through out the paper). 
It is only cut off, in an ideal world, by the macroscopic, or integral scale L. The probability density function 
(PDF) for the order parameter m associated with the diverging correlation length is essentially the exponential of 
the free energy P(m) ~ exp(— F(m)/k^T) and takes on an approximately Gaussian form as long as the Landau 
approximation, F{m) ~ a + bm 2 + is valid. Close to the critical point the Landau approximation breaks down and 
the PDF becomes non-Gaussian. The key assumption of the renormalisation group theory of critical phenomena is 
that the critical PDF remains scale invariant in the thermodynamic limit and can be obtained from the fixed point 
of a renormalisation group transformation J3|,|| (see App. ^). Thus, renormalisation group theory can be regarded 
as a generalisation of the central limit theorem to systems that are correlated over all length scales. The critical 
PDFs can be termed "universal", in that, when properly normalised, they depend on at most a few basic symmetries 
that define the universality class of the system. A non-Gaussian and universal PDF is therefore a direct signature of 
the fluctuation driven critical phenomena that have revolutionised modern statistical mechanics || . Analytical and 
numerical work on Ising, Potts, and XY models has shown that a generic feature of such systems is a skewness, 
with large fluctuations below the mean, towards small order parameter values. 

Correlations that are both strong and long range are a feature not only of critical phenomena but also of systems 
driven far from equilibrium. However, in the case of driven systems, the absence of a microscopic theory means that 
one has to rely heavily on empirical observations from experiment and numerical simulation. Labbe et al. ||To[| have 
shown that the PDF for the energy injected into a closed turbulent flow at constant Reynolds number is also non- 
Gaussian and universal. In this case "universal" means that the PDF, when suitably normalised, does not depend on 
Reynolds number or several other parameters (for example, the type of fluid). The PDF again has a marked skewness 
with an apparent exponential tail for fluctuations towards low energies. 

The present work is motivated by our empirical observation pd|Jl2] | that the universal PDF of energy fluctuations in 



the turbulence experiment 1 1 13|] is, within experimental error, of the same functional form as that of the universal 
P(m) for the critical system that we have studied. The latter is the spin wave limit to the low temperature phase of the 
2D-XY model [O 15) that is known to capture the critical behaviour of the full 2D-XY model Jl6| po[] . The distribution 
is shown in Fig7]2: it is asymmetric, with fluctuations below the mean approaching an exponential asymptote, while 
those above the mean approach a double exponential. This observation led us to the proposition that many systems, 
both equilibrium and non-equilibrium, sharing the property of long range correlations and multi-scale fluctuations, 
should share the same features, at least to a good approximation pi| , pi| . The proposition appears to be strikingly 
confirmed in ref. [pl|| where, from numerical simulation, similar behaviour is observed in a number of different systems: 
for order parameter fluctuations in the two-dimensional Ising model and in the two-dimensional percolation problem, 
as well as as for fluctuations in global quantities for models of forest fires and avalanches, driven into a self-organised 
critical state. This appears to contradict the idea that the PDF should depend on the particular universality class 
of the model in hand. One possible way of accounting for our observations is that many universality classes share 
common features, with the differences between them appearing either outside the range of physical observation, or 
being hidden by experimental error. There are therefore many open questions regarding a possible and much desired 
connection between critical phenomena and non-equilibrium systems as well as regarding the details of the PDF in 
critical systems. It is these questions that we address in the current paper, via an analytic study of the PDF for order 
parameter fluctuations in a finite XY model in arbitrary dimension. 



B. Normalisation of the Order Parameter 



We discuss order parameter fluctuations of finite systems in terms of distributions that are calculated in the 
thermodynamic limit, N — > oo. As discussed in App. |A|, it is essential to normalise the order parameter by an 
appropriate power of N = L D in order to obtain a distribution of finite width, or, equivalently, a form for P{m) that 
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is independent of system size. By extending the scaling hypothesis to include finite systems j22j, the following form 
of P(m) has been proposed: 

P(m, L) ~ L f3 / lJ P L {mL f3 / 1 ' ', f/L). (1) 

Here f3 and v are the conventionally defined critical exponents for the magnetisation and correlation length £ respec- 
tively p2[ . The appropriate normalisation of the order parameter is provided by the factor L~®l v while fixing different 
ratios £/L will in principle result in an infinity of different limit distributions as L — > oo. We concentrate on the case 
of a truly critical system with correlations over all length scales, which should result in maximum deviation from the 
Gaussian form. Here the dependence on £ can be dropped from eqn. ([!]), and Pl{iti) should closely approximate a 
single universal function of the variable mL^l v for all values of L. In this form it is independent of the microscopic 
details of the system, although it could indeed depend on the universality class of the transition through the critical 
exponents. 

Equation ([[[) is demystified somewhat by recognising that the normalising factor L~^l v is, in such an ideal system, 
proportional to the mean value of the order parameter, (m). Further, one of the key properties of a critical system is 
that the standard deviation of the distribution, cr, scales with system size in the same way as the mean value. This 
property, which is a direct result of the hyperscaling relation and which we refer to as the hyperscaling condition, 
means that eqn. ([!]) can alternatively be written in the form P(m) = l/aPi,{m/a). Thus a provides, as might be 
intuitively expected, the correct normalisation of the order parameter, such that a reasonable PDF of finite width is 
obtained in the limit N — > oo. In this paper, in addition, we shift the distribution with respect to the mean value and 
define 



aP(m) = 11(6) , (2) 

where 

= I^M. (3) 



In this representation one expects the PDF to fall, in the thermodynamic limit, onto a single universal curve. Provided 
that finite-size corrections to scaling are negligible one should observe data collapse onto this universal curve for large 
but finite system sizes. 



C. The Two Dimensional XY Model 



The model that we study, the harmonic spin wave limit to the XY model, is defined in section || for the case of two 
dimensions, D = 2. This is the dimension of most interest in the present context, as the system is at its lower critical 
dimension. At low temperature the coupling J /k^T is an exactly marginal variable that characterises a line of critical 
points in zero applied field [^3| . The critical line is separated from the paramagnetic phase by the Kosterlitz-Thouless- 
Berezinskii phase transition at Tktb ]l6| , |l7| ]. The critical phase that exists below this temperature is an attractive 
subject of investigation from both an analytic and a numerical point of view. Its physics is entirely captured by the 
harmonic Hamiltonian |l6|,[l9],^0| with the result that many calculations can be performed microscopically, without 
the need to use renormalisation techniques, or the scaling hypothesis. From a numerical point of view simulation 
results near a single, isolated, critical point are often complicated by a shift in the effective critical temperature by 
an amount scaling to zero as L^ x l v [0-|9), making it unclear exactly which temperature should be studied. Indeed 
numerical studies of Ising and Potts models |^-^| do find distributions whose form depends on temperature in the 
critical region (see App. [A]). In the 2D-XY model, as the system is critical over a range of temperatures there are no 
such technical problems and data for P(m) can be collected at all points below Tktb- These factors make the 2D-XY 
model an ideal system with which to study the effects of critical correlations. 

The finite-size scaling for the 2D-XY model has been discussed in our previous publications JlJ,|l5| . In this work, 
we began, following Berezinskii (l6) and Racz and Plischke (24), an exact calculation for the PDF of order parameter 
fluctuations. This calculation is completed and presented in detail in the current paper (section ||). It shows explicitly 
that the non-Gaussian behaviour in the 2Z3-XY model stems from the influence of all length scales from the microscopic 
to the macroscopic scale. We propose that the same is true for other complex systems including those driven far from 
equilibrium, which provides a basis for understanding the apparent overlap of their PDFs and provides an unexpected 
experimental motivation for studying a system as simple as the 2Z3-XY model. 

Two results coming out of our calculation are worthy of note at this stage. The first is an exact analytic result 
that is rather surprising given the previous discussion and the general belief concerning the dependence of the PDF 
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on universality class: shifting the curve with respect to the mean, eqn. (f2|), gives us universal data collapse, not 
only for all system size but also for all temperatures for which the harmonic Hamiltonian is valid. The ratio of 
exponents (3/v depends linearly on temperature, from which we deduce that the PDF is independent of the value 
of the exponents along the line of critical points. One should note however that these points are rather special and 
the result cannot necessarily be generalised to all critical points: not all the usual critical exponents are defined. For 
example the exponents (3 and v are not individually defined, but their ratio is |]l7|| . The usual scaling relations are 
valid in terms of the ratio (3/v only and this "weak scaling" |l8| means that there is only one independent critical 
exponent, r\ — 2(3 jv |rij| , compared with two for a regular critical point. This is all that is required for the analysis 
leading to eqn. , but is not sufficient to ensure a unique functional form for the general problem with two exponents. 
However, it does seem consistent with the idea that only small differences separate results for different universality 
classes. 

The second result that it is relevant to mention at this stage concerns the finite size scaling data collapse of eqn. (||). 
We find that the hyperscaling property a ~ (m) is not a necessary condition for data collapse onto a non-Gaussian 
function. With our definition (||), the first two moments of P{m) fall trivially out of the calculation and all that is 
required for data collapse is that the moments (9 P ) for p > 2, are independent of the system size. This is the most 
general condition for non-Gaussian data collapse, while the PDF only satisfies the scaling hypothesis in the form of 
eqn. ([!]), if the hyperscaling condition is satisfied. We give, in section II. A, an explicit example where data collapse 
onto the universal curve of the 2D-XY model occurs, but where the hyperscaling condition is not satisfied. If we 
make an expansion of the order parameter about a perfectly ordered state (m = 1) in powers of temperature, keeping 
only the linear term, then (m) diverges logarithmically with system size [ p5||2lj ], while the standard deviation is a 
constant. The ratio a/ (m) is actually an increasing function of system size throughout the physical domain. It is only 
when the order parameter is correctly defined on the interval [0, 1] that the hyperscaling relation is re-established, but 
written in the form (Q) the two distinct variables have the same universal PDF, even outside the range of temperature 
and system size for which the development gives an accurate representation of the true order parameter. This result 
is more than a mathematical curiosity; the harmonic approximation for the 2D-XY model maps directly onto the 
Edwards- Wilkinson (EW) model |^7|- |30| , ^2| for interface growth and the linearised order parameter is related to the 
square of the interface width, w: m — 1 — w 2 . Our PDF therefore corresponds precisely to that for interface width 
fluctuations and for which, in two dimensions, the hyperscaling condition for the observable w 2 is explicitly violated. 



D. Organisation of the Paper 



The rest of the paper is organised as follows. In section |[| we present detail of the calculation for the PDF in the 
2£)-XY model We show explicitly that it is a univ ersal function of system size and of temperature and find an 
exact expression for the characteristic function (section II A| ) . Transforming the distribution numerically, we compare 
it in detail with extensive Monte Carlo and molecular dynamics simulations of the full XY model and show that it 
is clearly the complete solution of the problem (section II B). We calculate the asymptotic values of the distribution 
for large devi ation s belo w and above the mean, which we find to be exponential and double exponential respectively 



(sections HC, and [ID) 



In section 



[II we try to fit the computed PDF to standard functions by comparing the moment expansion of the 



generating function with those of the Fourier transform of the test function. Three functions are considered: 



11(6) 



' exp a[9 — s — exp(0 — s)] , 

exp j- [log(s - 6>) - a} 2 j , 



(4) 



\v/2-l -a(s-e) 



The PDF is fitted to an excellent approximation over the physical range by the first two functions, while the thir d give a 
reasonable but slightly inferior fit. The first function, with a an integer, comes from extremal statistics (section [II A). 



For convenience, throughout this paper we use the term "XY model" to refer to either the model defined by the spin wave 
Hamiltonian or the full XY model over a temperature range in which the spin wave approximation is valid. This should not 
cause any confusion in reading the present paper, but our choice of terminology should be born in mind when comparing to 
other work on the XY model 
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It is Cumbers third asymptote, corresponding to the PDF for the a th largest value from ensembles of N random 
numbers |j33f . The interpretation with a non-integer (we find a = tt/2) i s not clear, but a connection between critical 
phenomena and extremal st atisti cs is a very appealing concept |34|,|35|]. The second function is a generalised log- 
normal distribution (section [II B ). Unlike the first curve, it does not have the correct asymptotic forms but despite 
this it fits just as well over the physical domain. The th ird function is a y 2 distribution describing identical and 
statistically independent degrees of freedom (section [ill C ) . It gives reasonable qualitative agreement indicating that 
a good, zeroth order description of a correlated system is in terms of a reduced number of statistically independent 
variables. However, this description has its limits, as shown by the fact that this function fits the exact PDF slightly 
less well than the other two. This variety of different fits suggests that one should treat the physical interpretations 
that they offer with caution; however even with this caveat in mind they still represent useful mathematical tools. 
To investigate this point further, in section [II D we derive an approximate functional form for the curve using an 
analysis due to Pearson which reconstructs the PDF from the four principle moments, which in this case have been 
calculated analytically. The Pearson analysis gives a quite different function which also gives a good description of the 
exact PDF over a physical range of 6. This serves to emphasise that, given zero mean and unit variance, the shape 
of the curve over a typical experimental range is essentially defined by its skewness, 7, and kurtosis, k. Therefore, an 
alternative way of summarising the observed universality [jn],^), is that 7 and k, for several different systems, have 
the same scale-invariant values as they do for the XY model. 

In section IV we extend our calculation to D dimensions, which apart from D — 2 are all non-critical. Despite this, 
we find evidence of the integral scale for all dimensions D < 4. For D = 1, the PDF for the linearised order parameter 
shows an exponential tail. However, we show numerically that the PDF for the correctly defined order p aram eter is 
quite different and is just what one would expect for a paramagnetic system without correlations (section [V A). The 
case D = 3 holds a final surprise (section [VB): despite the long range order of the low temperature phase the PDF is 
still not a Gaussian function. The temperature is a dangerously irrelevant variable in the ordered phase of the 3D-XY 
model with the result that the susceptibility remains weakly divergent at low temperature |36|. The result of this 
divergence is that the asymptotes of the PDF for large fluctuations are exponential below the mean and exp(— Cst 8 3 ) 
above the mean. The hyperscaling relation, in this case, is again violated. The divergence disappears at the upper 
critical dimension and we find a truly Gaussian PDF for D > 4. 

In section [v| we conclude by returning to the physical reasons for the exponential tail in the PDF. The XY model 
is diagonisable in reciprocal space reducing it to a model of statistically independent degrees of freedom: spin wave 
amplitudes at wave vector q, </> q . The amplitudes diverge at small q and are these modes that give the non- 
Gaussian fluctuations. In one-dimension they completely destroy magnetic order, in two-dimensions they give critical 
behaviour and between two and four dimensions they give remnant critical behaviour in the form of a dangerously 
irrelevant variable. 



II. PROBABILITY DENSITY FUNCTION FOR THE ORDER PARAMETER IN THE 2D-XY MODEL 



A. Analytic Expression 

The 2D-XY model is defined by the Hamiltonian 

H =-J^cos(^-6>j), (5) 

(i,3) 

where the angles di refer to the orientation of classical spins Si confined in a plane and where the sum is over 
nearest neighbour spins. In the following we consider a square lattice of side L, with periodic boundaries. The 
order parameter is a two dimensional vector m which, in zero field is free to point in any direction. We define the 
instantaneous magnetisation as the scalar m = |m| 

I N 

i=l 

where 9 = tan _1 (^ i sin^/ J2i cos#j) is the instantaneous magnetisation direction. Within small corrections, which 
disappear in the thermodynamic limit, this corresponds to the more conventional definition 
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For all temperatures below Tktb the renormalisation group trajectories flow, at large length scale, towards a regime 
where only spin- wave excitations are relevant 1 19 2(J . The physics of the low temperature phase is therefore completely 
captured by the quadratic Hamiltonian 



H 



(7) 



We therefore restrict ourselves, in the following calculation to this Hamiltonian and neglect the periodicity of the 
variables 9i. Our calculation can not therefore take into account the presence of vortex pairs. Close to Tktb in two 
dimensions and also in one dimension, where free vortices are relevant variables we would expect a deviation from the 
behaviour shown in Fig. [2| This point is discussed further below. 

We now calculate the PDF P(m) that the system be in a state with magnetisation m, using the standard property 
that a probability density function is entirely defined by the value of its moments j3?j]. Indeed, P(m) can be expressed 
in terms of its characteristic function, P(x): 



r°° dx ■ 
P(m) = I ^ e m *P(x) , 



(8) 



which can in turn be expanded in a Taylor series whose coefficients are the moments 



p(x) = E 



x p d p P 



p=0 



p\ dx p 



E 

p=0 



(—ix) 1 



(9) 



so that 



P(m) 



dx 
2^ 



OO 

; E 

p=0 



{-ixf 



(10) 



Eqn. (|T^) assumes that the series converges and that all the moments exists. Note that this last feature demands that 
P{m) falls off faster than any power-law of m. 

The program for calculating P(m) is therefore to calculate the moments (m p ), sum the series and transform the 
final result. To this end it is useful to define the Green function in Fourier space, 



G(q) 



1 



4 — 2 cos q x — 2 cos q y ' 



(11) 



where q x and q y take the discrete values ^-n of the Brillouin zone with n = 0,...L— 1. We also define the set 
of constants gk — Y^q G(q) k /N k . The value of g\ diverges logarithmically with system size, illustrating the critical 
nature of the low temperature phase [^,^6): g\ = ^rlog(CiV), C = 1.8456 |$8j. The values of gk, fc > 2 are 
independent of N in the thermodynamic limit. We find: g2 — 3.8667 10 -3 , (73 ~ 7.5719 10 -5 , 54 = 1.7626 10 -6 and 
that for large fc, gk behaves like (27r) 1_2,c /2(fc — 1), see App. 

The first moment is easily calculated within this approximation (see App. [b| and refs. [^[ll]]). One finds that (m) 
decreases algebraically with the size, as one would expect from finite-size scaling p2| 



(Ncy 



-k B T/8ivJ 



(12) 



As discussed above, while the critical exponents (3 and v are not individually defined for the 2D-XY model, their ratio 
is jl^l and the system obeys what Kosterlitz referres to as weak scaling [^8) . Through eqn. ( |l2"| ) the ratio of exponents 
is defined: j3/v = rj/2 = T/4nJ. 

For higher moments we need a more systematic approach. A specific property of the quadratic Hamiltonian (fj]) is 
that the moments can be calculated using the tools of Gaussian integration |l4|,^]. In particular, by the application of 
Wick's theorem, propagators of order 2p in reciprocal space can be exactly expressed in terms of quadratic propagators 
so that the p th moment is proportional to (m) p . One finds (l5| 



(2N)p 



E E 



exp 



=±1 



(13) 
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where r is the reduced temperature k-gT / J and Gr{t) the regularised Green function J2 q M ^(l) ex P(*Q 1 r )/N. In 
order to compute each moment of order p, we have to evaluate the sums over the positions and operators c^. The 
idea is to expand the exponential term ( |l3| ) and introduce a diagrammatic representation of each quantity computed. 
For example, we represent aiGu{vi — Yj)<7j by a line between i and j on a lattice of p sites. The general term of the 
expansion is then a set of graphs with a combinatorial factor for the symmetries. Since a 2 = 1, only closed diagrams 
are relevant, the factor TP being cancelled by the sum over all the 0$, The factor of r fe , is common to all graphs with 
fc lines connected together, with an even connectivity at each vertex. For example, up to the second order term in r, 
we have 



(14) 



The term ^ r G 2 R {r)/N 2 = Sq^o G(q) 2 /N 2 = g 2 is the value of the one loop graph with two lines, as shown in 
Fig. |a. There is an additional symmetry factor 2 x p(p — 1), which is the number of possible positions for such 
diagrams connecting 2 lines on a closed graph on a lattice of p points. For the third order term in t, we have only 
one diagram with 3 vertices, of value 33, Fig. 0b. The symmetry factor is equal to p(p— I)(p-2)x4x2. The factor 
4x2 comes from the number of possible ways of connecting 3 lines together. For the A th order term there are three 
different graphs, two of which are shown in Fig. 0c. The first has 3 loops and 2 vertices, the second, of value 174, has 
one loop and 4 vertices. The third graph, not shown, consists of two disconnected one loop graphs of the type shown 
in Fig. 0a. In general, at each order in r, we have the product of different closed diagrams, with one or many loops. 
It appears that the values of multiple loop graphs, like the first one in Fig. 0c), are zero in the thermodynamic limit. 
We therefore find that only the one loop diagrams are relevant and the value for such a diagram, with k lines and fc 
vertices is gk- We can now express the p th moment of the magnetisation as 



(m p } 



£ 

k>2 



I-! £ £ 



k\ 



9ki ■ ■■9kr C (ki 



r>l fciH \-k r = k 

ki>2 



, k r ) x p{p - 1) • • • (p - k + 1) . 



(15) 



with C(fci, . . . , k r ) a combinatorial factor which takes into account the possible ways of putting together k lines on 
r graphs, the first with k\ lines, the second with ki lines, etc., including the symmetries. For example, the factor 
associated with one triangle is C(3) = 4 x 2. It is then relatively easy to show that 



C(fci,. 



~ik—r 



k) 



(fci H h k r )(k 2 H h k r ) ■ ■ -K 



(16) 



Next, we can use the fact that every diagram is invariant by the action of the group S r of permutations of its r single 
elements, so that, instead of (161), one can use a more convenient form for the combinatorial factor: 



4 £ C{k a{1) ,...,k, 



a(r)) 



1 2 fc - r fc! 
r\ k\ ■ ■ ■ k r 



(17) 



Setting fk i — g^i (~ T ) ki /2fe» we arrive at the result 



(to) 



3 = i+EEti E fk 1 ---fk r xp(p-i)---(p-k + i) 



fe>2r=2 feiH yk r =k 

k 



exp 



fc>2r=2 fciH \-k r =k 

00 



,fe=2 



= 1 +EEi E fin-hrltf* 



(18) 



For p = 2 we find (m 2 ) / (m) 2 = 1 + g2T 2 /2 and defining a — y/ (to 2 ) — (to) 2 we thus arrive at the hyperscaling 
condition that the ratio a/(m) is independent of the system size. Hence 




(19) 
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One can now substitute for (to p ) in (|h]) using ([TJ]) and after re-arranging the summations the distribution can finally 
be expressed as an integral, depending on the values of the one loop diagrams gk only 



P(m) 



dx 
2^ 



exp 



ix (m — (to)) + — ^(ir(m) 



fe=2 



Changing variables, x — > x/a and using jis| ) we find 



f°° dx 



TO — (m) 



ST 9k f ■ /2 



fc=2 



(20) 



which is the principal result of ref. [TBI. Defining aP(m) = n(#), we see that the function II depends uniquely on 
the variable 9 = (to — (m))/a and the gk, k > 2. As the ^ are constants in the thermodynamic limit, 11(0) is a 
universal function, independent of both system size and temperature. The asymmetry comes from the fact that the 
ratios 3fc/(<72/2) fe ^ 2 , k > 3 are non zero and this constitutes the abnormal influence of the integral scale. If, in the 
thermodynamic limit, k = 2 were the only non-zero term one would arrive at a Gaussian PDF centered on (to). 
Departure from a Gaussian function is typically characterised by the skewness, 7 = (9 3 ) and kurtosis, k = (9 4 ) [(4l[ . 
We find 
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93 



(ff 2 /2) 3 / 2 

34 



3 + 3 



(ff 2 /2) 2 



-0.8907, 
4.415. 



(21) 



a) 




b) 



"A 



c) 




q+k+1 




FIG. 1. Diagrams contributing to the distribution: a) to order r 2 , b) to order r 3 and c) to order r 4 . 

Although we can calculate the asymptotic behaviour of gk for large fc, we are not able to compute the constants 
analytically and so we cannot sum the series (pp|). However we can transform it into a very much more useful form by 
keeping N large but finite and inverting the sums over q and k. The even and odd terms are separated and summed 
independently and we eventually find: 
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11(9) = 



Fglidx i . iy 2 
— — — — exp < ixu \ I — — > 
2 2tt v 1 V 2 ^ 

q#0 



~zG(q)/iV - i arctan (xG(q)/N) + I log (l + x 2 G(<i) 2 /N 2 ) 



(22) 



The sum over q and the integral over x, in (E2J) can now be performed numerically, allowing the evaluation of H(0). 



B. Comparison with Simulation 

To test the above calculation and to verify its scaling properties we have carried out extensive numerical simulations 
of the 2D-XY model with full cosine interaction, eqn. (|J), for different values of temperature and system size. In 
addition, we have also done microcanonical molecular dynamics (MD) simulations to check the possible dependence 
of the PDF for fluctuations on the statistical ensemble. 

The Monte Carlo simulations were performed with 10 s Monte Carlo steps per spin, with 10 6 steps used for equili- 



bration. The MD simulation was carried out for systems of N classical rotators, 42 , with Hamiltonian 

N 2 

^ = Ey + J E[ 1 - cos ^-^ ■ ( 23 ) 

i=l (ij) 

The equations of motion were integrated numerically, using a Verlet algorithm. In order to explore the low-temperature 
fluctuations regime the initial configuration of the system was chosen with the spins pointing in the same direction 
and with a Gaussian distribution of momenta. The system was then equilibrated for a time of 10 6 — 10 7 sweeps 
and data collected over a time span of 10 8 — 10 9 sweeps according the size of the system. Note that one cannot use 
the harmonic interaction ^ to study deterministic dynamics in the microcanonical ensemble, as this would allow no 
coupling between the spin wave modes and no evolution would be possible. The non-linearity of the cosine interaction 
allows mixing between the normal modes and the sampling of equilibrium states. Here we do not report work at high 



enough energies to allow vortex formation 43 44J] with any significant probability. Rather, the non-linearity plays the 
role of the heat bath in the canonical ensemble, while the physics is still correctly described by the harmonic part of 
the interaction. 

The numerical integration of eqn. J22]), performed with a fast Fourier transform (FFT) algorithm fill , is shown in 
Fig. |^, where it is compared with Monte Carlo results for T/J = 0.1 and N = 32 2 . The theoretical curve is clearly in 
extremely good agreement with the numerical data. The curve is asymmetric, with what appears to be an exponential 
tail for fluctuations below the mean, with a much more rapid fall off in amplitude, for fluctuations above the mean. 

In Fig. U we show the PDF for fluctuations in m obtained from MC simulation for fixed system size and varying 
temperature, as well as MD for fixed temperature and different system sizes. The result of ref. Jll|] and section [n] of 
this paper is that, for the harmonic Hamiltonian, eqn (R), 11(0) is independent of both system size and temperature, 
while we have explicitly tested this result against the PDF generated for the full Hamiltonian, eqn. (||). Qualitative 
agreement is clearly excellent, independently of the ensemble used, but there are small systematic deviations in the 
tails, when observed on a logarithmic scales [[l6|, as shown in Figs. |i| and |5|. We can only expect agreement between 
the analytic result and simulation in the range of temperature sufficiently below Tktb such that vortex pairs do not 
influence the PDF fl4|| . Even in the absence of vortices one must expect small variations from our theoretical result 
for small system sizes that stem from the utilisation of (^|) rather than (pi). In a renormalisation group treatment 
the non-linearities of Hamiltonian (|^) scale away on changing the length scale and the Hamiltonian is replaced by 
an effective harmonic Hamiltonian at higher temperature [ pp[ . For example, at T/J = 0.7, for L = 32, we find 
(rn) = 0.76 from simulation, while eqn. (|12j) gives (rn) = 0.81. The effective coupling constant can be calculated by 



expanding the cosine and approximating the nonlinear terms using a Hartree approximation 19 1. Renormalisation 
of the non-linearities introduces a microscopic length scale a' which gives small corrections when compared with the 
calculated PDF. However, this length scale is fixed by the temperature and the corrections should scale away as the 
ratio a' / L — > 0. This scenario is confirmed in Fig. ||, where data are shown at T/J — 0.7, for L = 8, 16, 32, 64 and 
compared with the theoretical curve. Deviations from the theoretical result are observed for L = 8 and L = 16 but 
the PDF clearly approaches the predicted scale independence for the larger system sizes. 
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FIG. 2. The PDF, as obtained from a fast Fourier transform (FFT) of equation (^), compared with MC simulation of a 
2D-XY model at temperature T = 0.1 of size iV = 32 2 (upper: natural scale; lower: semi-log scale). 

Near Tktb vortices influence the PDF, however the vortex population decreases exponentially moving away from 
^ktb p3| ] and they only make their presence felt within the physical domain in a small band of temperatures near 
the transition. In this regime the data do not fit on the universal curve [|lj,[44|,[l(| but a detailed discussion of this 
point is outside the scope of this paper. 
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FIG. 3. The PDF for fluctuations in dimension D = 2 from MC and MD simulations. The first set of data corresponds to 
canonical MC simulation for a system of size N = 32 2 at temperature T = 0.1, 0.3, 0.5. The second set of data corresponds to 
microcanonical MD simulation at temperature T ~ 0.7 and size N = 16 2 , 32 2 , 64 2 . 
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FIG. 4. The PDF, as obtained from a fast Fourier transform (FFT) of equation (|22j), in dimension D — 2 compared with 
Monte Carlo results for a system of size N = 32 2 at temperature T = 0.1, 0.3, 0.5, 0.7. 
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FIG. 5. The PDF, as obtained from a fast Fourier transform (FFT) of equation (p2|), in dimension D — 2 compared with 
MD results for a system of size N = 8 2 , 16 2 , 32 2 , 64 2 , at temperature T ~ 0.7. 



C. P(m) for the Linearised Order Parameter 



As eqn. (|22j ) is independent of temperature, one should be able to obtain it at low temperature where the magneti- 
sation is approximately 



27V ^ 



(ft - of 



(24) 



In fact, using this expression one can arrive at (|22|) in a more straightforward manner. What is perhaps surprising is 
that the calculation, using (|24|) is valid for all temperatures below Tktb, even for temperatures where (^) and ( |24] ) 
represent different physical quantities. 
Using the Hamiltonian (Q), we have 



where is the inverse Green's function operator connecting sites i and j with non-zero elements for i and j nearest 
neighbours Jig ], and Z = (det G _1 /27rr) is the partition function. 

It is easy to integrate the Gaussian integral by transforming into reciprocal space. Defining the trace Tr of any 
function of G as the sum for q ^ of the same function of G(q) and and using (to) = 1 — tTtG/2N we find 



P(m) 



dx 

2^ 6XP 



T 1 

ix (to - (to)) - ix-TrG/N- -Trlog(l - ixrG/N) 



(25) 



We can now use the fact that a = \J gi /2t in this approximation, to transform ( p5j ) into a dimensionless and universal 
form 



11(0) 



— OO 

oo 



fg^dx 
2 2tt 

2 2tt 



exp 
exp [i$(a: 



/.r« 4 / — - i^TvG/N - ^Tr log (1 - ixG/N) 



(26) 



which is the same expression as (|2C],|2^), once we separate the real and imaginary parts of the integrand. This 
demonstration proves that the only relevant graphs are those with only one loop, the others being zero in the 
thermodynamic limit. 
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Within this linear approximation the mean magnetisation (to) and the standard deviation a do not scale in the 
same way with system size: while (to) = 1 — (T/8irJ) \og(CN)j_a — y / <72/2r is a temperature dependent constant. 
This exact result can be verified by applying eqn. (^|) to eqn. ( p5| ) and calculating (to) and (to 2 ) directly. The fact 
that we find the same universal function for the two calculations, when written in the form (0) shows explicitly that 
the hyperscaling result, a / (m) ~ 0(1) is not a necessary condition for non-Gaussian data collapse. Rather, it seems 
that hyperscaling is a consequence, in these circumstances, of the correct definition of m as an order parameter on 
the interval [0, 1]. 

The Gaussian limit of the 2D-XY model is identical to the Edwards- Wilkinson model of interface growth and the 
linear approximation for the ord er p arameter is related to the square of the interface width to = 1 — w 2 . The PDF 
for w 2 has been studied in one p8fl and two pi| dimensions together with extensions to the EW model, including 
non-linearity All models give non-Gaussian PDF's with the same qualitative features as Fig. ||. These models 

provide an important microscopic link between equilibrium and non-equilibrium systems and suggest that a formalism 
could exist that incorporates the statistical features that we have observed to be shared, at a global level, between 
such different systems. 



D. Asymptotes of 11(9) for Large Fluctuations 



As a first step towards an analytic form for 11(0) one can approximate (|22j) beyond the Gaussian approximation by 
retaining only the elements (<72,<?3)- In this case, the solution is proportional to the Airy function 



11(0) oc exp I - I Ai 
6a 



1 



(3o)V3i2a (3a)V3 



(27) 



where a = 2 3 / 2 g 3 /3gl /2 ~ 0.296876. The g 3 term assures that it is not symmetric on reversing the sign of 9. We find 
that the approximation reproduces qualitatively the apparent exponential behaviour for <C — 1: 



11(0) ~ < 2^ 



-,1/4- 



(3a)V3l2a (3a) 1 / 3 



exp ■ 



1 

6a 



1 



(3a)V3l2a (3a) 1 / 3 



3/2' 



(28) 



However the approximation does not allow us to extract the asymptote above the mean, as for 9 > the Airy function 
develops oscillations. 

A more fruitful approach is to look at the saddle points of the integrand ((2^), from which one can extract both 
asymptotes. If 9 <C —1, an expansion near x = is not very satisfactory and we must rather seek the solution for the 
extrema of the whole integrand, d<fr(x)/dx = 0. We find: 



^=I T r G 



iV 3 1 + x 2 G 2 /N 2 



1 G 2 x 

2 N 2 l + x 2 G 2 /N 2 ' 



(29) 



If 9 is negative and x real, the real part of the second term is always positive and there is no solution to this equation. 
We therefore seek a solution for x pure complex, x — iy. In this case, eqn. ( p9|) becomes 



—9 = -Tr ^ 



N 2 1 + yG/N 



(30) 



The function ip has simple poles at y = —4ir 2 , Sir 2 , —32ir 2 , . . . and its asymptotic value near the first pole yo = —Att 2 
is f(y) ~ —2/(y — yo). The extremum of the integrand satisfies the condition y* ~ yo — 2y/2/g2/9 > yo, for \9\ large 
and we can deform the real path of the integration so that it passes through the extremum on the imaginary axis. 
Near the extremum, we can expand the integrand up to second order in y — y* and perform a Gaussian integration: 



11(9) 



32 dx 
2 2~^ 



exp 



.1 



i$(iy*) + i-(x - iy*) $ (iy*) 



We finally find that the asymptotic value of the distribution varies as 

11(0) cx |0| exp ( 4n 2 ' 92 ' 



(31) 



(32) 
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We have superimposed the asymptotic result ( |32] ) and the full numerical integration for N — 101 2 of (22) in Fig. ||. 
The amplitude of eqn. ([32]) is chosen so that the curves are slightly displaced, to allow comparison of the slopes. 
The asymptotic solution is in excellent agreement even for 9 values where the PDF shows a distinct deviation away 
from exponential behaviour and only fails for 8 > —2. Further out in the tail, in the range — 10 < 8 < —4, log(II) is 
approximately linear. However the value of the slope is not the argument of the exponential in (j32|), 47r 2 -y/(?2/2 ~ 1.736. 
The logarithmic corrections given by the the term \8\ are significant over the whole of this range, but the curvature 
is so small that the data can be fitted to an effective exponential H(8) ~ exp(a#), with a — 1.56867. . .. The data 
only approaches true exponential behaviour for 8 < —30, which is completely outside any imaginable physical range. 
Strictly speaking it is therefore more correct to speak of pseudo-exponential, x exp(ax) , for the asymptote below the 
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FIG. 6. Comparison of the tail of the PDF with the exact asymptote (long dashed), eqn. (^2|), the true exponential tail of 
slope 4-7T 2 g2 jl — 1.736 (dotted), and an effective exponential tail of slope a = 1.56867. . . (short dashed). The curves are 
displaced from each other for clarity. 



For large and positive 8 a solution of eqn. (|30|) exists for large and positive y. A reasonable approximation is to 
replace G by 1/q 2 and perform the integration 

1 Nd 2 q 1 y 

~ 2 J q=27T/VW 4,r» 7VV 1 + y/Nq* (33) 
dq 1 

7T lo sv 



4tt J 2 n/^v e(i + <i 2 ) 8?r 

A more precise computation gives ip = log(y)/87r + a + l/2y + ■ ■ ■ where o is a numerical constant which can be 
computed exactly. An analytical study (see App. |c|) gives 

1 1 1 °° 

a = 24 + 4 " 47 l0g(4?r) " 27 l0g II t 1 cx P(-2^)] - -0.11351444337 . . . (34) 

fe=i 

For large 8, the saddle point of the integrand is therefore located at y* = exp87r(— a + yf^fff), and the asymptotic 
value for n follows from a Gaussian integration of (pl|): 



11(8) oc exp 



8tt V 2 



(35) 



Comparing the asymptote with the full curve we again find that the true asymptote only fits accurately outside the 
physical domain, although the data is clearly consistent with a very rapid fall off in the PDF for 8 above the mean. 
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III. FITTING TO KNOWN FUNCTIONAL FORMS 



The obvious question now arises: is the PDF generated by the characteristic function (E2|) of known functional 
form? We do not have a definitive answer to this question, as we are not able to transform (|22|) analytically. In the 
absence of an answer, we test the PDF against three skewed functions, shown in eqn. (jij), which describe statistics 
in different physical situations. These are: a modified Gumbel function, characteristic of problems where extreme 
values dominate the sum over many contributions; a log-normal distribution, characteristic of statistically independent 
multiplicative processes and a x 2 distribution which describes the PDF of a quantity made up of a finite number of 
positive definite microscopic variables. The analysis is the same in all three cases, but is only shown in detail for the 
modified Gumbel function: each curve has 4 parameters, but once the value of the first is chosen the others are fixed 
by normalisation and the constraints (9) = 0, (0 2 ) = 1. The family of one parameter curves are Fourier transformed 
and the first four terms in a Taylor expansion are set equal to those for the generating function, which fixes the 
value of the free parameter. The method takes into account the skewness of the curve but not the kurtosis and its 
accuracy is ultimately limited. The goodness of fit can be measured by comparing the ratio of higher order terms of 
the expansion of the test and generating functions. For an exact solution all higher ratios would be equal to unity, 
while for a poor fit they diverge rapidly from this value. Other functions could be tested in the same way and an 
exact solution may well exist in the statistics literature, unknown to us. 

The method described above is quite similar to that due to Pearson 49 1, who realised a century ago that, in all 
practical situations, knowledge of the first four moments of a distribution is sufficient to generate a curve, fitting any 
set of data points [p2| . Pearson developed a phenomenological differential equation containing the numerical values 
of the moments, whose solution gives the fitting function. A Pearson analysis is performed on the calculated PDF at 
the end of the section. 



A. The Generalised Gumbel Distribution 



The asymptotes (|32|) and p5| ) are of the same general form as those for Gumbel's first asymptote distribution from 
the theory of extremal statistics J33|: defining z to be the a th largest value from a set of Zi, i = 1,N random numbers 
taken from a generator f(z); the PDF for z is 



9a{z) 



^ exp {-a [a a (z - u a ) + e"^-^)] } . 



r(o) 



(36) 



r(a) is the gamma function; u a is the value of z such that a of the N random numbers are greater than z. F(z) is 
the probability of having a of the values less than z, such that F(u a ) = 1 — a/N. a a is referred to as the intensity: 
a a = (N I a) f (u a ) . In conventional statistics, a would of course be an integer. However, in what follows we are going 
to see an irrational number appearing. 

The function (|3^) has an exponential tail for fluctuations towards large values of z, the opposite of the PDF, in 
Fig. |^. We therefore make a change of variables m z = 1 — z, Z = (m z — (m z ))/a z which makes a mirror reflection of 
eqn. (|3^). Within the linear approximation for the order parameter this corresponds to the relevant variable being 
the sum of the spin wave amplitudes z — > (1/2 AT) J2i(@i — ?_1 5lj]. Changing variables we find 



v t n G (9,)=we<*°'-')-'>* K ''-* 
b = a a a z 

s = (1 - (m z ) - u a )/a z 
a a a a 



w 



r(a) 



(37) 



Eqn. ( |37j ) is also the distribution for the a th smallest random number from the set z^. After some algebra one can 
show that 




<9 2 r(a) 



da 2 



1 dT{a) 



log(a) 



T(a) 
0T(a) 



da 



r(a) da 



(38) 



Now re- writing (32) and (pq), one finds 
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|0|exp (§60) , 9 < 

11(6) oc { ' ' (39) 

exp [-^e b{ - 9 -^ +c6] , 6 > 0, 



with & = 8ny/g 2 /2 — 1.105, s = 0.745 and c — b. These asymptotes differ oniy slightly from those for a generalised 
Gumbel function with a = 7r/2: firstly through the term 16*1 for fluctuations below the mean and secondly through 
the term exp(c0) above the mean: the coefficient c = (%/2)b for the modified Gumbel function, while c = b for the 
true asymptote. These differences are enough to ensure that the modified Gumbel's equation is not an exact solution 
to eqn. (pp|), however the comparison is so close that it is tempting try to get a good fit to ( p2|) by solving for the 
constants a, b, s and w. 

Fourier transforming (k57]) gives 



dx iv / x \ / x \ f dx 

— — exp ( ixO — isx + i— loga — aloga I T I a — i— ) = / — exp [i&c(%)] ■ (40) 
-00 2lT ° v ° J \ bJ J_ 00 2tt 



We can compare <&g(x) with $>(xy / '2/g 2 ) 1 assuming that the two Fourier transforms are nearly equal. The four 
constants should be calculated by minimising the difference between the two functions. To do this we can set the first 
four coefficients of the Taylor expansion of these functions equal. For $>g(x) we have: 

$g(x) = ialoga - i\og(w/b) - ilogT(a) + [s - *(a)/& + log(a)/&] x + -^S> '(a)x 2 

where ^f(z) is the digamma function V (z) /T(z). For <& we have: 

*W27£) = \^ - ^ - + + (42) 

L 3g 2 < A9 2 hg 2 ' 

We therefore find that the four constants satisfy the relations: 

b T(a) 

- = ^T' S 6 = loga-*(a), 

& 2 = *'( a ), b3 93 ^Y = -*» (43) 

The first three equations arise from the constraints of normalisation of the distribution, while the last expresses these 
constraints in terms of g 2 and (73. The equations can be solved numerically. We find 

a = 1.5806801, 6 = 0.9339355 

s = 0.3731792, w = 2.1602858 (44) 

The constants b and s calculated in this way are shifted slightly from the values extracted from the asymptotes, 
but a is close to our very appealing first try ir/2. Taking this value and calculating the constants &, s and w from 
normalisation one finds: 

a = tt/2, b = 0.938 

s = 0.374, w = 2.14, (45) 

in very satisfying agreement with the first method of calculation. 

Given this solution, we can compute the coefficient ratio for the higher order terms in ( |4T| ) and d42|): 



1 d 4 <t>(x^2/f 2 ) 



$W(0) a* 4 



g^"'(a 

x — 



= 0.9265029, (46) 



1 d 5 ${x^/2/gz) 



$(5)( ) dx 5 



48V2 55 6 5 _ 0826742g (4?) 



gl ,2 ^)(a) 
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The ratio of coefficients clearly does diverges from unity, but is does so slowly, indicating that the modified Gumbel 
function should be a good fit to the curve over the physical range. This is confirmed in Fig. [7] where we compare (|37| ) , 
using the values with the exact result, from eqn. (^||). On a natural scale the agreement is remarkably good over 
the entire range, with the only visible deviation coming around the maximum of the PDF, where the Gumbel curve 
is very slightly lower. On a logarithmic scale there is excellent general agreement over the whole of the plotted range, 
but a slight deviation can be observed for probabilities below 10~ 3 . For fluctuations below the mean the deviation 
is because the true asymptotic behaviour is quasi-exponential, x exp(— ax) and has a slight curvature, as discussed 
in the previous section. The results therefore confirm that, although the generalised Gumbel function is an excellent 
approximation for the PDF (|22|), it is not an exact solution. 

^,From these results it is very tempting to take the generalised Gumbel function, with a exactly 7r/2 as a working 
analytic expression for the PDF. However the connection with extremal statistics remains an open question p5| . 
As discussed in section [v], the spin wave Hamiltonian (Q) is diagonalised in reciprocal space and the problem can 
be formulated in terms of a set of statistically independent variables. The PDF for extreme values of statistically 



independent variables can only follow three different asymptotic [|33 41|, or limit functions as the thermodynamic 
limit is taken. The only possible limit functions from extremal statistics of the Gumbel form discussed here are for a 
integer, with a = 1 for the biggest or smallest values. 

Chapman et al. have recently argued that the PDF for global quantities in any system with identifiable 

excitations on scales up to the system size should be dominated by extreme values. They shown that the PDF of 
extreme values among 10 Gaussian random number generators approximates to a Gumbel function with a = tt/2. 
This is not one of the predicted asymptotes EIL and wc suggest that the deviation must be due to a very slow 
approach to the limit function with system size. It therefore does not seem to be a correct description of the 2D-XY 
data as we do have a limit function which is well represented by eqn. ([57]) with a — tt/2. However, if the results 
of |35| ] are relevant for non-equilibrium phenomena such as turbulence and self-organised criticality it would suggest 
the interesting property that corrections to the asymptotic forms, or limit functions, are a generic feature of these 
systems. 

B. Generalised log-normal Distribution 

The generalised log-normal distribution has the form 



with w = 1. Following the same procedure as before, the generating function <&l(x) can be developed as a power 
series 



(x) =x(e-S + c Q +^/ 2 ) + i^-(e 2a + 2 °t ~ e 2a +^) - x 3 (l e 3a+9*l/2 + l_ e 3a+3*l/2 _ i e 3a+5^/2^ _ (4Q) 



LW = 

Comparing (|4^) with (Q) one finds the following expressions for s, a and <tl 

„a+£r|/2 



S = C 

a = -iln(e 2,T ' - e CT ' ) 
= I e 3a fe 9 ^/ 2 + 2e 3CT '/ 2 - Ze 5 °H 2 ) . (50) 

3 g^ 1 6 V / 

Eliminating a and ol leads to a cubic equation for s in terms of a = g ^ — = 1/ItI : 

s 3 - 3as 2 -a = 0, (51) 

which could be solved exactly. We have solved it numerically, verifying that there exists one real and two complex 
roots. We find 

s = 3.45981, a = 1.20109, a L = 0.28325. 
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The function, with these parameters, is compared with the calculated PDF in Fig. |7[ The general quality of fit is 
again excellent over the plotted range, with very small systematic deviations occurring in the wings of the distribution. 
It does not have the correct asymptotes; either exponential on the left or double exponential on the right, but as we 
have shown in the previous section, the true asymptotic behaviour is only reached outside the plotted regime, which 
explains why such a good fit can be achieved. 

We have not, for the moment been able to develop a physical reasoning associated with the log- normal function 
and the origin, s = 3.4... although related to 7, seems rather arbitrary, but we do not exclude an explanation in terms 
of random multiplicative processes. 

Note that log-normal distribution does appear in surface dynamics. Namely, starting with a flat interface as an 
initial condition, the short-time limit of the D — 1 Edwards- Wilkinson dynamics yields a log-normal distribution for 
the interface width p0[ . 

C. Generalised \ 2 Distribution 

The x 2 distribution for v statistically independent degrees of freedom has the form 

U x (6) = w(s - oy/ 2 - 1 *-*'- 8 ), (52) 

with 

a"' 2 



I>/2) 

v = 2a 2 . (53) 
As in the case of the Gumbel function, the generating function can be found in closed form: 

$ x {x) = x{6 - s) +i-ln(l -ix/a), (54) 

whose development up to 4 th order in x leads to 

* x (x) = x{6 - s ) + ^x + i^x 2 - £^x 3 - i-^x 4 + 0(x 5 ) + ... (55) 

This series can again can be compared with ( |42| ) to give 

v 

s = — 
2a 



2 9 • 
S3 T 



(56) 



with numerical values 

v = 10.07155, a = 2.24405 , s = a ,w = 2.31233. 

Comparing the function, shown in Fig. [?] with these parameters with the calculated curve, there is reasonably good 
agreement but this time deviation can be seen when plotted both on real and logarithmic scale. On the logarithmic 
scale the deviation is stronger than for the other fitting functions. 

One can see that describing the correlated system as a finite number of degrees of freedom is a reasonably good 
approximation. It is an appealing concept and the calculation yields a system size independent number which depends 
uniquely on the skewness: v = g\jg\ — 8/7 2 . If 7 developed towards zero, then v would diverge and the \ 2 
interpretation would be consistent with a Gaussian distribution. However, quantitatively it is not correct and the 
true description is a many body one (5^] . The difference between the two curves can be quantified by considering the 
ratio of the 4 th order terms: 

= -il 

* (aj )W = -<|L > (57) 
so that ^(2) (4) ~ 0.0238 which is very far from 1. 
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FIG. 7. The PDF compared with the generalised Gumbel, log-normal and \ 2 functions described in the text. 

D. Pearson's Curve 

Pearson jl9],[30| described an ingenious method of deriving a functional form for a PDF to fit experimental data, 
given the first four moments of the latter. He considered the differential equation 

— ; — = ; (58) 

ax bo + b\x + oar 

and showed that if y is a distribution then the parameters b, bo, b\ are specific functions of the first four moments. 
The expression can then be integrated to give (within a normalisation factor) an approximate functional form for the 
PDF, which by definition has the same principal moments as the data to be fitted. The success of Pearson's approach 
relies on the observation that PDFs with the same moments are approximately coincident over the range of a few 
standard deviations, which is exactly the range of experimental interest. In the present case the mean is zero and the 
standard deviation is set to unity, so the shape of the curve depends only on the skewness, 7, and kurtosis, k. 
We find 7 = g 3 (2/g 2 ) 3/2 = -0.8907, and n = 3 + 3st 4 (2/g 2 ) 2 = 4.415, which gives the following solution: 
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(59) 



in which £ = x - 0.39723, (5 = 2.4787, a = 11.430, q = 10.249, p = 47.267, y = exp (105.02). Equating y(x) = U{6), 
the fit to the exact expression (Fig. |8|) is good between x — — 6 and x = 2, but the very large numbers involved in 
eqn. ( |59| ) suggest that this functional form has no physical significance. From this analysis we can conclude that data 
collapse observed in refs. ]2l}| should be interpreted as meaning that the third and fourth moments scale with a and 
L in the same way as they do in the critical 2D-XY model. 
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FIG. 8. The PDF compared with the fit obtained with the Pearson method described in the text. 

IV. DISTRIBUTION IN THE D-DIMENSIONAL GAUSSIAN MODEL 

In this section, we study the asymptotics of the distribution function in general dimension D. It is straightforward to 
generalise the development from eqn. ( po|) to eqn. ( |22] ) for arbitrary dimension by redefining G(q) for dimension D and 
summing over a D dimensional Brillouin zone. The generalised expression (|22|) can then be numerically transformed 
to give II((9). The results for D = 1 and D = 3 are shown in Figs. [l(] and [11|, where they are compared with data 
from Monte Carlo and molecular dynamics simulations. There is again excellent agreement showing that eqn. (|2^) is 
essentially exact in the low temperature regime, where the Hamiltonian (]?]) is valid. At higher temperatures the full 
Hamiltonian (||) generates vortex structures, the eqn. ( |l3| ) is no longer valid and the derivation of eqn. ( |2^ ) breaks 
down. Within the low temperature approximation there are three regimes: D < 2, 2 < D < 4, and D > 4, in addition 
to the special case D = 2. The different regimes can be seen from a dimensional analysis of g\ and gi. As deviation 
from Gaussian behaviour is due to the abnormal influence of the integral scale in the form of infrared divergences, 
we approximate replacing G by 1/q 2 and re-calculate the gu by performing integrals over the Brillouin zone between 
2-k/N 1 / d and 2it. This procedure gives the correct qualitative behaviour, but there is a difference between the discrete 
sums and the integrals over the Brillouin zone, even in the thermodynamic limit (see App. |^). The correct qualitative 
behaviour is 

{c hD M 2 -°y D , D < 2 
AilogiV + Bi, D = 2 (60) 
C x ,d , D>2 

and 



20 



r c 2 , D N 2 ^- D y D , d < a 

52 -< (A 2 logN + B 2 )/N, D = A (61) 
I C 2 ,d/N , D > 4. 

The lower and upper critical dimensions, D = 2 and 13 = 4, are marked by the logarithmic divergence of g\ and 52 
respectively. 

Using the linearised order parameter (|^) we find for D < 2 that <?i diverges as a power of N giving (m) = 
1 — tCi i dN( 2 ~ d ^ d , which is a poor approximation for a thermodynamic quantity bounded on the interval [0, 1]. 
Once outside this restricted low temperature region, r < 1/ [Ci^AT^ 2- - ^ 1 , both the linear approximation for the 
order parameter and the quadratic approximation for the Hamiltonian break down and there is a divergence in the 
behaviour of the PDF, as calculated from (j^) and as simulated numerically. The system is, of course, disordered at all 
temperatures, so that the correct (to) and a both vary as 1/yN and the PDF for the vector order parameter is a two 
dimensional Gaussian function centered on m = 0. The PDF for to, as defined in (||), is P(m) ~ mexp (— m 2 /2a 2 ), 
analogous to a Maxwellian distribution of molecular speeds, and the thermodynamic system satisfies the central limit 
theorem (see App. |a|). As we have already seen, for D = 2 the situation is different, as there is a large region of 
temperature where the quadratic Hamiltonian correctly describes the physics even though eqn. (0) is not a good 
approximation. In this regime of temperature, the PDF 11(0), for parameters (§3) and (Q) are however identical. 

For dimension D > 2, the low temperature expansion for the order parameter gives consistent results for all N, as 
long range order is stable and (to) ~ 1 — Ci^r is well defined. Above D = 4, our results agree with mean field theory 
(D = 00) where all sites are connected together. Here, (m) ~ 1 — r/4 and a ~ r/2v / 2A", and for large but finite N, 
the universal function n is simply a Gaussian 

n(*) = ^exp(-^), (62) 

which corresponds to the central l imit theorem for a collection of N independent oscillators, each of expectation value 
(to) and standard deviation r/2y2N. 




FIG. 9. The PDF in one dimension (N = 128) at temperature T/J > 12/N. The continuous line is Maxwell speeds 
distribution of an ideal gas. 
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FIG. 10. The PDF in one dimension (N = 128) at temperature T/J < 12/JV". The dashed line (with slope ~ 1.04) is the 
exponential asymptote for the low temperature approximation given by eqn. (p5j) and is shifted with respect to the main curve 
for clarity. 



A. Low Temperature Calculation in D = 1 



If the low temperature calculation for D < 2 is not terribly pertinent for the thermodynamic system, it it highly 
relevant for the interface problem in the context of the EW model |^,^8|,^9) and is exactly solvable in D = 1 |2^] . In 
this case, computing the different gu, we find 



, 91 = TV/12, g 2 = A 2 /720, . 



2C(2 P ) 



p > 1 



where £(fc) = * fe is the Riemann zeta function The expectation value of the magnetisation and standard 

deviation are 



(to) = exp 



tN 
"24" 



1 - 



tN 
~2T' 



coshrA(x — x + l/6)dx — 1 I (to) 



t z N z 



which means that the ratio (m)/er scales as 1/A, although for the parameters of the interface problem 
(w 2 )/a w 2 ~ O(l). We evaluate the universal distribution II, using the general eqn. ( p6| ) with G defined for the 
D = 1. After some algebra, we find for H(0) 



n / gexp 

dx 



2vr 



exp 



/360 

/360" 
12 



$>g 1 



k=l 



log 



zxv360 

27T 2 fc 2 



sin W ixv360/2 



/ 2^ GXP 



^^360/2 



(63) 



This expression is related directly to the function $ of eqn. (11) in ref. pq| : 11(0) = $(2 — 246>/\/360). The method 
used in |^4|,^8],^9| is based on path integration, but the results are the same as our saddle point method, used to 
compute the asymptotics. Setting x = iy, the extrema of $ satisfy the equation 
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^/y-v/360/27!-2 



C0th7T 




For 9 <C — 1, y is close to the first pole — 27r 2 /V360 of the right hand side of Q), which is similar to the 2D case ( |30| ) 
except that the ID extrema function is easier to evaluate. Performing the saddle point computation, we find that fl 
behaves asymptotically as 



11(6) oc exp 27r 2 6>/\/360 



(65) 



which is the same as p8j. The asymptote, (65), is drawn on Fig. (|10| b), where it can be compared with the full 
calculation and with simulation. The exponential tail is extremely well defined and the predicted slope is clearly 
correct. 

In the regime of fluctuations above the mean, for y ^> 1, 6 is close to the constant V360/12, and no extrema exist 
for 9 beyond this value. In this case, y ~ V360/[8(V360/12 — 9) 2 ], and the saddle point approximation leads to the 
following asymptotic value for II near this upper limit 



11(0) cx 




9) 



(66) 



. We refer the reader to ref.s 



for the precise coefficients in both asymptotic 



which is the same result as 
limits. 

In conclusion, we find that for 9 -C — 1, the universal distribution again has an exponential tail, while for large 
fluctuations above the mean the PDF shoots to zero as exp[— 39q/2(9 — 9 )], with 6>o = V360/12. This upper limit 
corresponds to the constraint that m < 1. 

It is worth pointing out in some detail here that the exponential tail in the one-dimensional problem is not the 
result of critical fluctuations. The small deviations in angle (9i — 9j) constitute a random walk with w ~ \/l — m 
being the radius of gyration, which scales correctly as the square root of the walk length, L. The ID linear order 
parameter, or interface problem is therefore nothing more than a simple random walk p8j, but despite this the PDF 
is, as shown in Fig. |l^: a standard result for such a walk is that the mean radius of gyration is proportional to the 
mean end to end distance, S of the walk. It is easily shown that the PDF P(S) is Gaussian Changing variable 
from S to X — S 2 one finds P(X) ~ X -1 / 2 exp(— X/Xq); a trivial distribution with an exponential tail. The PDF for 
w 2 has the same exponential tail, but does not show the essential singularity at w 2 — (m — 1) and we conclude the 
rather surprising property of a random walk, that the PDF for the radius of gyration and for the end to end distance 
are not the same. The origin of this difference is that the average angle 9, corresponding to the center of mass of 
an equivalent random walk, fluctuates with L in the same way as the radius of gyration itself and this lack of self 
averaging removes the essential singularity from the PDF at w 2 = 0. 



B. Asymptotic Solutions in General Dimension 

We first evaluate the asymptotic value of II for positive 9 by solving the saddle point of (|2^) , rescaling the variable 
X V 52/2 — > x for convenience. For D < 2, the ratio 91 1 ^[92 is independent of the system size and, with x = iy, the 
equation to solve is 



.91 



Cst 



Nq 



D-l 



Cst /jvi/d q ^N + y^f2jg: 



zdq 



where Cst is a constant. By setting N 1 ' D q/ ^Jy — > q, we find that, for large and positive y 



9i 



cx y 



(D-2)/2 



Cst N^ D /^y D-l 



Cst /y/y 



l + q' 



dq^y 



(D-2)/2 



00 q - 1 
l+q 2 



dq 



(67) 



(68) 



which means that 9 is close to the upper bound g\j '^/2g~2~. Replacing the asymptotic value of y for the extrema in the 
function <f> (l2r3), we find that 
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, D/(D-2) 

9i 



logll(0) Cst -p=-0 + Logarithm corrections, 6 p=, D < 2. (69) 

V V 2 92 / V 2 ff2 

The logarithmic corrections come partly from the Gaussian integration around the saddle point and partly from other 



terms in (68) which are not accurately evaluated within our approximation. Note again that D = 2 is a special case 
as, instead of (|6^) we have a logarithmic divergence (see eqn. (p3|)) and subsequently a double exponential fall in II 
for large 9. For the interval 2 < D < 4 the ratio g\j ^ftyi and the integral ( |68| ) are no longer finite and so we look to 
eqn. (J30|) for the asymptotic behaviour: 

1 f Cst d D q y 



oc 



92 Jest /NV° Nq 4 1 + yVZ/i^Nq 2 ) ' 



(70) 



By again setting N X I D qj ^Jy — ► g and using the fact that giN 2 ^ D 2 ^ D is finite (60), we arrive at 



■/o 1 + T 

The integral is convergent for 2 < D < 4 and by replacing the value for y in the saddle point approximation, we get 
the asymptotic form for II, in the limit of large and positive 9: 

log 11(0) Cst 9 D/[D - 2) + Logarithm corrections, > 1, 2 < D < 4. (72) 

In 3 dimensions, we therefore expect that the logarithm of the distribution falls off like 9 3 , well above the mean. We 
have not tested this in detail, but the PDF does fall off more slowly for D = 3 than D = 2, in qualitative agreement 
with the predictions here. Finally, we note that throughout the range 2 < D < 4 the universal PDF is non-Gaussian, 
but the hyperscaling relation is invalid: (m)/a ~ giZ-i/gi ~ N( D ~ 2 ^ D . 

For D > 4, gi decreases as 1/N, consequently, (JTl|) has to be modified. We find, instead of (fn]), that 

9 (X / £ „ yj N >> L (73) 

J A r(D-4)/4D/yy 1 + g 

We can, in fact, replace the integrand inside the integral by q D ~ 5 dq since the integration domain is large, from which 
we find that the saddle point is proportional to 9 ^> 1 and deduce that LI is Gaussian on the right hand side of the 
curve. The same is true for D — 4 despite the logarithmic divergence of 52- 

In the opposite limit 9 <C —1, for both D = 1 and D = 2 the asymptotic value of the distribution falls down 
exponentially ( |32] , |65|) . We would now like to evaluate this limit in general dimensions. In both cases the coefficient of 
9 is related to the value of 32, i.e. 6*2,15. Rewriting the eqn. ( p0| ) with discrete sums (see also App. |c]), we have 

Q _ \ a f 74) 

16 ^ ff2 ^0 K + -"+™d) K + ... +TO 2 3 )+y^27^iV(2-C)/^/4 7r 2 

where the sum excludes rrij = 0, i = 1, ■ ■ ■ , D. The saddle point equation has a solution y which is the pole nearest 
the origin, y — — 4n 2 ^ g2/2N( D ~ 2 ^ D , i.e. for sets of {mi} with one element equal to 1, the others being zero. For 
D < 4 and large N, this pole is finite since (72 compensates N 2 ( D ~ 2 ^ D , so that its value is simply y = —4tt 2 ^/C2~d/2. 
Applying the saddle point integration, we find that the dominant term in the logarithm of II is, below the mean 



logIT(0)~47ry^0, £<4, (75) 

and is linear in 9 for every dimension below 4. Included in Fig. [n] for D = 3 is a fit, on the left hand side of the 
form (75), with (^2,3 calculated numerically. There is again excellent agreement, which convincingly confirms the 
presence of the exponential tail. In fact, true exponential behaviour is reached for smaller values of 9 than for D — 2. 

For D > 4, the value of this pole diverges like jV( D ~ 4 )/ 2D , and the previous solution fails. In fact, the solution ( [73] ) 
for positive 9 and y is also valid for negative values, if q D ~ 1 dq/ (q 2 + I) is replaced by q D ~ 1 dq/(q 2 — 1). Since the 
integration domain is far from the pole of the denominator, we can approximate the integrand in both cases by 
q D ~ 5 dq, and we get the same result as (]7^). We therefore, finally conclude that LI is also Gaussian on the left hand 
side of the curve and the central limit theorem applies for D > 4. 
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FIG. 11. The PDF in three dimensions for N 



asymptote given by eqn. 



and T/J = 1.82. The dashed line (with slope ~ 2.5) is the exponential 



and is shifted with respect to the main curve for clarity. 



V. CONCLUSION 



Probability functions with exponential, rather than Gaussian behaviour are a common feature of complex sys- 
tems |54],fl7 48 5^-p6|. For example, the PDF for velocity differences at microscopic scales, in fully developed turbu- 
lence show exponential tails J54|. This appears to be true in turbulence, not only for microscopic quantities, but also 
for global quantities; the energy injected into a closed turbulent flow being a very well controlled and documented 
example M , ^3 57 1. Following these observations we have proposed that this is also a generic feature of complex sys- 
tems [ p"l|]21| ]. In this paper we have shown that, for the low temperature phase of the XY model, a critical system at 
equilibrium, analogous behaviour occurs when a few long wavelength and large amplitude modes make their presence 
felt in the global measure, which is typically a sum over O(N) degrees of freedom. The exponential tail can occur in 
three physically different situations. The first is in two-dimensions, when the system is critical and fluctuations occur 
over all length scales. The second is in one-dimension, when the system is not critical, but an exponential tail occurs 
for a particular global measure, relevant to problems of interface growth, whose moments are completely dominated 
by the integral scale. The third is in three-dimensions, also non-critical, where despite stable long range order, the 
large amplitude long wavelength modes continue to make their presence felt. The detailed form of the PDF in these 
three cases are quite different and easily discernible in experiment. In table 1 we show the evolution of the skewness 
and the kurtosis with spatial dimension. The deviations from the Gaussian limit are largest in one dimension, and 
decrease continually to zero at D = 4. We propose that the difference in the form of the PDF could be used as an 
experimental signature of the underlying physics. 

From the general evolution shown in table 1, one might expect a dependence on shape, with dimensional crossover 
as the length scale in one direction changes from microscopic to macroscopic. This is indeed the case and for example, 
in two-dimensions, the skewness and kurtosis of the PDF calculated from eqn. ( p2| ) increases towards the values for 
D = 1 if the ratio of lengths in the x and y directions, L x and L y , are varied continuously from unity. It would be 
extremely interesting to establish if the same is true when the length scales are varied in turbulence experiments and 
numerically, in the models of self-organised criticality. 

To see how the anisotropy of the PDF comes from the long wavelength excitations, we give an analysis in reciprocal 
space: the Hamiltonian (0) is diagonalised 



(76) 



q>0 



where 4>q is the discrete Fourier transform of 8i 
variable for each q taken as the real part of <j), 
written m = 1 — Y]^ 



and the sum is over the Brillouin zone [p8|, with the thermodynamic 
,. Defining m q = (l/2iY)5ft{0 q } 2 the linear order parameter can be 



iq, where the m q are statistically independent variables with PDF 
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P ( TOq ) = yW^ TOq l/2 e -/3^ mq _ (77) 

Here, as we are principally interested in the modes at small q — |q| we have, without loss of generality, approximated 
G(q) _1 ~ q 2 . The PDF for to is thus nothing more than the composite PDF for a set of independent spin wave modes 
or an "ideal gas" of particles, whose only peculiarity is that the mass term varies as q~ 2 . The Goldstone modes have 
wave vector q = 2n/L and hence make contributions of 0(1) to m, while the modes on the zone edge, with q = 7r have 
only microscopic amplitude. This dispersion in amplitudes is the key to the unusual behaviour for D — 2, as it violates 
one of the conditions for the central limit theorem to apply to a sum of statistically independent variables: that the 
individual amplitudes do not differ by too much. However it is not true that the Goldstone modes, by themselves 

give the complete PDF. The mean value (j2q m q^ ~ I^-k/lI 271 ^^ wnere ~ q D ~ l is the density of states. 
For D — 2 both limits of the integral are required and a detailed calculation gives (^2 q m q ) = (ry/4) log(CiV), with 
C = 1.87 and with critical exponent r\ = T/2nJ. The anomalous term logiV therefore reflects the fact that modes 
from all over the Brillouin zone are relevant for (to) and through eqn. ([l8]), for the higher moments (m p ). 

For D = 1 only the lower limit of integration is required, the upper limit can be set to oo and the constants g p 
are proportional to N p . As a result the linear development of the order parameter in small angles, (Pj), is a very 
poor approximation for the thermodynamic quantity defined on the interval [0, 1]. The two expressions,"^]) and (|24|), 
describe different physical quantities. The former is directly related to the interface width in the Edwards- Wilkinson 
model of interface growth. The PDF for the full order parameter is consistent with an uncorrelated system, that is, a 
paramagnet with two-dimensional order parameter. For the linear order parameter the PDF, shown in Fig. |l^, does 
have an exponential tail, but this is not the result of critical fluctuations, it is the property of a simple random walk. We 
remark further that dependence on a macroscopic length scale does not, in itself, indicate critical behaviour. Rather, 
critical behaviour is exemplified by the case D = 2, where all length scale are important between the microscopic and 
macroscopic cut off. 

D = 3 represents the opposite of the one dimensional case: (to) is controlled by the upper limit of integration and 
the result is unchanged by setting the lower limit to zero. However, despite long range order being stable and the 
system not being critical in the low temperature phase, the exponential tail persists. This is related to temperature 
being a dangerously irrelevant variable |59[| near the zero temperature fixed point of a renormalisation group flow, 
between the lower and the upper critical dimension. The constant g2 now falls to zero with system size but it does 
so more slowly than 1/N (see eqn. (|6l|)). As a result of this slow decay, the ratio g p /g^ 2 , p > 2 in eqn. (|2(]) is 
independent of N and the distribution is non-Gaussian, despite g p and gi both being zero in the thermodynamic 
limit. At low temperature the magnetisation is finite, but the Goldstone mode influences the PDF sufficiently to 
produce an exponential tail. A physical consequence of this anomaly is that the longitudinal susceptibility 

X~% ((m 2 ) - (to) 2 ) ~ N^' D (78) 

is weakly divergent throughout the ordered phase |l4|,|36|]. This is true for all magnetic systems with Heisenberg 
or XY symmetry. It could therefore be interesting to look for evidence of the departure from Gaussian behaviour 
experimentally in a non-critical three-dimensional system. Precision temperature control would not be required, 
however, as the ratio <t/(to) falls off as 1/7V 1 / 3 , the divergence in the susceptibility is very weak and this phenomenon 
may be out of experimental reach. 
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TABLE 1 Variation of skewness 7 and kurtosis k with dimension D 

Returning finally to critical systems; we have been able to exploit a system interacting via a quadratic Hamiltonian 
at exactly the lower critical dimension. In this particular case one has access to a critical point, with the fluctuation 
dominated behaviour that this implies, while retaining the benefit of Gaussian integration over phase space. As a 
result, all critical behaviour can be calculated microscopically, without the need for either the renormalisation group 
or the scaling hypothesis. The only price one pays for this simplicity is a critical system with a single independent 
exponent and the scaling relations satisfied through weak scaling only. In general, we believe that the analytic results 
that we have obtained are useful for the understanding of finite-size scaling and for the interpretation of experimental 
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observations from more complex correlated systems. The examples we have discussed Jll| , [2l| point towards a behaviour 
analogous to criticality for an enclosed turbulent flow and for models showing self-organised criticality. However the 
detailed analysis presented here leaves many open questions and more experiment and simulation are clearly required 
if the generality and the limits of this proposition are to be tested further. 
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APPENDIX A: 

Some comments on the central limit theorem in critical systems 

The central limit theorem is a powerful result of probability theory that provides the foundation for statistical 
thermodynamics |6(|. It states that the PDF of the sum Z = YliLi z i °f N statistically independent variates Zi 
tends, in the limit of large N and for moderate values of the variate Z ', to a Gaussian distribution. As well as 
the statistical independence of the z,-, another key criterion for the theorem to hold is that the Zi are individually 
62] ]. At a critical point, the first of these criteria is violated. The 2D XY model is of particular interest 
here as it is diagonalisable into statistically independent degrees of freedom and maps directly onto a problem where 
the second criterion is violated: the direct space variables, that is, the spins Si, are certainly individually negligible 
for large system size N, but are strongly correlated. On the other hand, when diagonalised in reciprocal space the 
spin wave variables are statistically independent, but are no longer all individually negligible. In particular, the long 
wavelength modes make a significant impact on the fluctuations of the global measure; in this case the linearised 
order parameter (|24|). The PDF for the full and the linear order parameters are identical, even when the quantities 
themselves differ, which makes it an ideal system for the practical study of the breakdown of the central limit theorem. 
A conventional critical system cannot, in general be diagonalised in this way, as evidenced by the divergent specific 
heat. 

Strictly speaking, the central limit theorem does not apply to the compound variate Z, but rather to the normalised 
quantity (Z — (Z^/N 1 / 2 . This normalisation is essential for a reasonable PDF in the thermodynamic limit, as the 
standard deviation for fluctuations about the mean value (Z) scales with system size in the same way. If a normalisation 
factor A fl / 2+p , p =/= 0, is chosen then one obtains a distribution that is concentrated either at zero or infinity S. We 
illustrate this with an example from statistical thermodynamics. The total energy E of an ideal gas of N molecules has 
a PDF of the form P(E) ~ I? 3Ar / 2_1 exp (—(3E). It is straightforward to confirm that P(E) tends to a delta function 
in the thermodynamic limit, while P(E /TV 1 / 2 ) tends to a Gaussian function pjfl . One can see from this example that 
the function is never truly Gaussian - indeed it is always of the form InP ~ (3N/2 — 1) In. 23 — 0E, which can easily 
be made independent of N by choosing appropriate units. The central limit theorem applies because the width of the 
distribution scales as N 1 / 2 which means that fluctuations with any physical significance are all concentrated near the 
turning point of the function In P. The theorem only has meaning because of the significance one attaches to values 
of the variate that differ by only a few standard deviations from the mean. In practical terms it is therefore essential 
to normalise fluctuations to the standard deviation in order to test the central limit theorem. 

In the case of dependent variables, the limit distribution can be different from the Gaussian form. Two types 
of dependent random variable can be defined |}J: (i) weakly dependent, in which the correlation function falls to a 
constant value in a finite range, and the standard deviation again varies as \/N ; (ii) strongly dependent, in which the 
fluctuations vary as a power of N different from one half. Case (i) corresponds to a system with a finite correlation 
length. In case (ii), which includes systems with critical fluctuations, the central limit theorem does not hold, but a 
reasonable PDF can be obtained by normalising to the variance, hence to an appropriate power of N, with p ^ 0. 
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Defining the (scalar) order parameter to be the intensive quantity z = Z/N, and using the scaling relations for a 
finite system, one finds p = (1 — rj/2)/D. The limit distribution is now expected to be non-Gaussian, as can be shown 
explicitly for the Ising model [^||34j. Note however that, p remains non-zero even at the upper critical dimension 
(taken as D = 4 here), when n — and where one might legitimately expect a Gaussian PDF. The condition p ^ 
may therefore be a necessary but not a sufficient condition to ensure non-Gaussian order parameter fluctuations. 

Case (ii) is not actually limited to critical fluctuations: the example of a dangerously irrelevant variable discussed 
in the text also falls into this category, with p = 2/D — 1/2. Here, p does go to zero as the upper critical dimension 
is reached and the danger of the irrelevant temperature variable disappears. 

An ordinary critical point is more complicated than those of the 21?- XY model. In this case the correlation length 
is only infinite precisely at the critical temperature. A non-Gaussian limit function can therefore only be found on 
a locus of points such that £/L is a constant as the thermodynamic limit is taken. Thus, fixing the temperature 
T =/= Tc and varying N will always cause a transition from non-Gaussian to Gaussian statistics. Conversely, fixing 
T = Tc one will only arrive at the stable limit function in the thermodynamic limit. One can therefore imagine 
a set of loci of constant PDF in [T, L^ 1 } space that converge on [Tc,0]. We have suggested |^) that there is one 
such locus, [T*(X),L _1 ], where the PDF has approximately the same form as that of the 2D-XY model. Thus, to 
sit at the critical temperature and change L is not the same as traveling along the locus [T*(L), L~ l \, From scaling 
argument |j5) one can check that the tails of the PDF at Tc should have the form P(m) ~ exp(— m s+1 ) in order to 
yield the correct scaling relation in the presence of a weak magnetic field: (m) ~ h 1 ^ . We do not find this, despite 
the same scaling relation holding for the 2D-XY model with S = 8irJ /ksT — 1. This difference may come from the 
difference in trajectories in the space of variable T and L. 

A final point concerns the central limit theorem as applied to a vector order parameter, m, such as the XY model. 
In the high temperature limit, the fluctuations in the vector m follow a two-dimension Gaussian centered on m = 
and the PDF for the scalar m = |m| follows a "Maxwell speed distribution" for a two-dimensional gas. In an ordered 
regime and even in the critical regime for D = 2 Q, a < (m) which means m behaves, to an excellent approximation, 
as a one-dimensional quantity. The symmetry breaking therefore induces a change in topology for the fluctuations in 
m. This is generalisable to order parameters of higher dimension. 

APPENDIX B: 

The graphs <?& can be written, in the large N limit in terms of power series. For example: 

, 4 1 4*1 
o 2 = km _ _ \ \ _ _ 1 _ \ (Bl) 

a^oo N 2 (4- 2cos2vrm/V]V - 2 cos 27m/ VA 7 ) 2 n ^ (4 - 2 cos 2irm/y/N) 2 



where Q — (v N — l)/2. The sum is dominated by the contributions for small m and n, but as the pole m = 0, n = 
is explicitly excluded from the sum, it remains finite even in the limit N — > oo. Taking only the first terms in a 
development of the cosines, which is exact in the thermodynamic limit one finds 
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4^2 2^ ^4 + ^4 ( m 2 +n 2)2 
m— 1 rn — 1 n— 1 / 

^ ^ oo oo 

m— 1 n— 1 y 



and in general, for 



-j^ -j^ OO OO -j^ 

gk ~~ 47T^ ^ m 2k + 47t4 ^ ^ (™ 2 4- rtW ' ^) 



m— 1 



(m 2 + n 2 Y 

-1 n— 1 v J 



APPENDIX C: 



For large and positive y the functional form of <p is ip(y) ~ logy + constant. To evaluate it in detail we use the 
results of the App. M to write: 
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V( y ) 0-2 E f ra 2 m 2±J" l "'uEE( m 2± n 2 ™2.„2 .,",)' ( C1 ) 



Q^oo 2tt 2 ^-^ \ m 2 m 2 + y J 2tt 2 ^— ' ^— ' V to 2 + n 2 to 2 + rt 2 + y 

m— 1 v 7 m— 1 n— 1 v 

where y = y/Air 2 The first two summations give, in the limit of large Q, a constant and a function of y which tends 
to zero for large argument: 

Hm J_ y ( i _ _J_^ = 1 (C2) 

y^oo 2lT 2 ^— ( V 777^ 771 + 7/ / 12 
m— 1 v ' 

The double sum can be rewritten as 

2 -2 > )=Arf f-^-?- 2 \ -R(Q,y) (C3) 
2^2 \ m 2 _|_ „2 m 2 + n 2 + yj 2?r 2 ^ ^ \ m 2 + n 2 m 2 + n 2 + y J v / 

m— 1 n— 1 v 7 m— 1 n— 1 v 7 

where i? is a correction term which vanishes in the limit of large Q: 

1 Q oo , 1 v 

fl(Q,») = rjE E 2 ^ 2X7 ( C4 ) 

277^ \ 777^ + 77^ 777/ + 77^ + y I 

ra=ln=Q+l v y/ 

The sum can be evaluated in the limit Q — > oo using the Abel-Plana formula |]66| : 

£ /(•) - jf /«* + i/W + i/M + 2 f ^'^y^ (C5, 

where / is any real function that satisfied the assumptions in ]6^ ]. Applying this to R{Q 7 y), we have 

1 1 .„ f°° xdx 1 



^-^4E^ W - -tan( Q /777)] + — + 4Q / 



m— 1 

i_ 

4Q ' 



2(m 2 + Q 2 ) (a; 2 - m 2 - Q 2 ) 2 + AQ 2 x 2 exp(27ra;) - 1 



^ = ir/2 — arctan (q/\J m 2 + y 

[ a / m 2 -4— L \ 

xdx 1 



1 



2(m 2 + y + Q 2 ) 



(C6) 



/ (a: 2 - 777 2 - y - Q 2 ) 2 + AQ 2 x 2 exp(27rx) - 1 
The first term tends, in the large Q limit, to the integral 

EQ 1 f"^ dx loff xdx 
— [it/2 — arctan(<5/777,)] — > / — [77/ 2 — arctan(l/x)l = / 5- = Catalan. 

m=1 TO Jo x io 1+^ 

A similar behaviour is found for the fourth term, since in this limit the dependence on y of this term vanishes as 
y/Q 2 . The other terms ar e cor rections proportional to the inverse of some power of Q, so that R vanishes in the large 
Q limit. The double sum (C3) can thus be reduced to a simple sum, since 

V „ 1 „ = --^ + COthTTZ. (C7) 

^ n 2 + z 2 2z 2 2z v ' 

n—l 

We therefore have, for large y 

1 1 1 \ 1 A 1 77 u 1 

7~77 > > 5 o o 5 7 = 77-77 > — 77 n + ~ COth 77777 - 

>tt 2 \ m 2 4- n 2 m 2 4- n 2 4- n / 9-n-2 ^-^ 9m 



277 2 ^— ' ^— ' V 777 2 + 7l 2 777 2 + 77 2 + V J 2tT 2 ' 2?77 2 2m 2(t77 2 + «) 

m-ln-1 v 7 m— 1 v x 



77 



2\J m 2 + 7/ 



COth 77 yV+^ 



77 + E "i T=r> ~ f 1 ~ coth77\/m 2 +y) + (coth 77777 - 1) 

24 ^— i 477A/777 2 4- 7/ V / 477777 

(C8) 



24 ^ 477 v^T^ 

±(L 1— 

477 \ 777 ,/m 2 +y 
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The series containing the hyperbolic function of y vanishes i n th e limit of large y and the asymptotic behaviour of 
the last term can be evaluated with the Abel-Plana formula (05) 



I rr 



i \ r j / 1 i \ i „ r xdx 

dx\ , + - + 2 



m — 1 



»' y/m 2 + y J Ji \x y/ x 2 + y I 2 J (1 + x 2 )(exp 2ttx - 1) 



= log (l + v / l+£) - log 2 + 7 (C9) 

where the constant 7 is equal to 



1 „ f°° xdx 
7= t: + 2 



2 J (1 + x 2 )(exp27rx- 1)' 

This can be proved by again applying the Abel-Plana formula to the function 1/m, since we know that Y^m=i V TO : 
logn + 7. The constant — coth7rm)/47rm in (|C8|) can be rewritten as 



£ 4^ (GOth " m ~ 1} = £ 2^ £ exp(-2 7 rmn) 

m— 1 m— 1 n— 1 

00 

= -^logII[ 1 - ex P(- 2 ™)] ( C1 °) 



and finally the results (C£ ,C£ , C1C ) give the asymptotic behaviour of ip for large y: 

<p(y) = i log j, + i - i- log 4^ + ^ - i- log f] [1 " exp(-2™)] + i- + • • ■ (Cll) 

n=l y 

The last term comes from a further study of the Abel-Plana formula which gives the other correction terms in the 
inverse power of y. An identical analysis gives the finite size magnetisation 



exp(-^TrG/7v) (C12) 



where TiG/N can be expand as 



i-TrG= -^logCW 
N 4tt 



C = exp <| I + 2 log — + 2 7 - 4 log [1 - exp(-27rn)] [ = 1.8456 (C13) 
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